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p . We examine the peculiar velocity distribution function of galaxies in cosmo- 

Q . logical many-body gravitational clustering. Our statistical mechanical approach 

derives a previous basic assumption and generalizes earlier results to galaxies 
c§ ! with haloes. Comparison with the observed peculiar velocity distributions indi- 

cates that individual massive galaxies are usually surrounded by their own haloes, 
^ ! rather than being embedded in common haloes. We then derive the density of 

energy states, giving the probability that a randomly chosen configuration of 
galaxies in space is bound and virialized. Gravitational clustering is very efficient. 
The results agree well with the observed probabilities for finding nearby groups 
containing N galaxies. A consequence is that our local relatively low mass group 
is quite typical, and the observed small departures from the local Hubble fiow 
beyond our group are highly probable. 
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1. INTRODUCTION 

A fundamental question about any spatial configuration of galaxies is whether it forms 
a gravitationally bound group or cluster. If so, we can ask whether it is sufficiently relaxed 
to satisfy the virial theorem, 2{K) + (W) = 0, relating the time averages of its kinetic and 
potential energies. A closely related problem is to estimate the effects of configurations on 
the peculiar velocities of galaxies around them. 

Approximate answers to these questions are known for computer simulations (Saslaw 
1987, 2000) and some observed clusters (Barychev et al 2001). In this paper, however, we 
are concerned with statistical properties of many configurations rather than individual cases. 
This will help determine whether the local group, say, is a typical or unusual configuration. 
Therefore we will explore probability distributions for peculiar velocities, total energies and 
virialization in regions of gravitational clustering. Since the observed clustering is highly 
non-linear on small scales, our theoretical description will also be fundamentally non-linear. 

The peculiar velocities of galaxies, i.e. their departures from the global cosmic expan- 
sion, provide a basis for our discussion of binding and virialization, and contain significant 
information about the history of galaxy clustering and the geography of dark matter. Much 
of this information can be represented by the velocity distribution function f{v)d^v, which 
is the probability that a galaxy has a peculiar velocity between v and v + dv. In the case 
of a perfect gas, this is the familiar Maxwell-Boltzmann distribution. In the case of galaxy 
clustering, which is highly non-linear on scales less than several megaparsecs, f{v) d^v de- 
parts greatly from the Maxwell-Boltzmann form. This is because it includes galaxies with 
all degrees of clustering, from isolated field galaxies to the densest rich clusters, in a very 
large representative volume of space. 

The present strongly non-linear velocity distribution presumably developed from more 
quiescent initial conditions. Initial states which start with a linear Gaussian distribution of 
density and velocity perturbations can be followed reasonably well into the weakly non-linear 
regime using perturbation theory (Nusser and Dekel 1993a,b; Bernardeau et al 2002). On 

scales where the distributions have evolved into strongly non-linear systems, dynamic dissi- 
pation has destroyed most of the detailed memory of the initial state. For the cosmological 
many-body case, this relaxation on small spatial scales, where the gravitational field is very 
grainy (Saslaw et al 1990), takes only one or two expansion timescales. These non-linearities 
then spread to larger scales as the clustering evolves. 

The observed velocity distribution function was discovered (Raychaudhury and Saslaw 
1996) using a representative sample of galaxies from the Matthewson survey which has rela- 
tively little a priori bias regarding their degree of clustering. Previously it had been predicted 
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(Saslaw et al 1990) using a quasi-equilibrium gravitational many-body theory of galaxy clus- 
tering. This theory in which galaxies may be surrounded by individual haloes, is in excellent 
agreement with relevant N-body simulations and agrees very well with the observations 
(Crane and Saslaw 1986; Fang and Zhou 1994; Shcth ct al 1994; Saslaw and Haquc-Copilah 
1998). Other models, such as those dominated by large cold dark matter haloes containing 
many galaxies could also in principle compute their resulting galaxy peculiar velocity dis- 
tribution functions and compare them with the observed distribution function, but they do 
not yet appear to have examined this question in detail. 

This distribution function in velocity (or momentum) space complements the distribu- 
tion in configuration space, and their consistent combination provides a statistical description 
of evolution in the complete six-dimensional phase space. Originally, (Saslaw et al 1990), 
the velocity distribution was derived from the spatial distribution by making the additional 
main assumption that on average the potential energy fluctuations are proportional to the 
kinetic energy fluctuations in any local volume. Direct A'^-body simulations confirmed this 
as a good approximation. Recently a new and more general statistical mechanical derivation 
(Ahmad, Saslaw and Bhat 2002) of the spatial distribution function has been calculated. 
In the present paper, we show how the earlier assumption relating fluctuations can now be 
proven directly from the statistical mechanical partition function, thus deriving the velocity 
distribution more rigorously. Moreover, the statistical mechanical derivation contains an 
explicit softening parameter for the gravitational potential. This provides a simple "isother- 
mal sphere" model for individual galaxy haloes and allows us to determine the effects of 
such haloes on the peculiar velocity distribution. Comparison with observations gives an 
approximate upper limit to the size of such haloes. 

Combining the velocity distribution function with the new statistical mechanical ap- 
proach also leads to relatively simple solutions of some other fundamental questions. The 
additional questions we will examine here are: What is the probability that N galaxies in a 
randomly placed volume of size V form a bound cluster? What is the probability that such 
a bound cluster is virialized? What is the probability of forming a small group of galaxies 
around which the peculiar velocity dispersion is so small that the local flow agrees well with 
the global Hubble expansion, as is observed around our local cluster? 

To answer these questions, we calculate the density of energy states for the gravita- 
tionally interacting cosmological many body system. This is one of a very few statistical 
mechanical interacting systems (the Ising model is another) for which such a calculation has 
been possible. 

In §2, we derive what was previously the basic assumption that local kinetic energy 
fluctuations are proportional to local potential energy fluctuation. We then obtain the ve- 
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locity distribution function for the softened potential representing simple isothermal haloes 
around individual galaxies. Comparison with the observed peculiar radial velocity distribu- 
tion function provides an estimate for the maximum average size of such haloes. §3 calculates 
the density of energy states and the probabilities for bound and virialized clusters, and §4 
discusses implications for the cool local Hubble flow. Finally, §5 briefly summarizes and 
discusses our results. 



2. DERIVATION OF THE VELOCITY DISTRIBUTION FUNCTION 

INCLUDING GALAXY HALOES 



To derive the velocity distribution, we start with the canonical partition function derived 
for the cosmological many body problem viewed as a model for galaxy clustering (Ahmad, 
Saslaw and Bhat 2002; Saslaw 2003): 



Zn{T, V) 



1 
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5^^ + '^('^i'^2,...riv) 



(2.1) 



where A is a normalization factor, to be discussed in §3, 



(2.2) 



assuming all N = nV galaxies in the volume V have an average mass m and temperature 
T, and 
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Here e is the softening parameter for the gravitational potential energy 

0(ri,r2,...riv) = ^ 0(rij-) , 

l<i<j<N 



with 



0j 



(2.4) 



(2.5) 



and Ri is a comparison scale which may, for example, be the average distance between those 
galaxies which contribute most of the mass. 
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When Vij < e, the potential energy is approximately constant, as would characterize an 
approximately isothermal sphere with density p oc r~^. Of course, this can be generalized 
easily with a range of values of e for galaxies of different masses, but here we are interested 
in basic results which can be obtained more simply by considering all the galaxies to have 
their average values of e and m. 

To examine fluctuations in the numbers and energies of galaxies which can move between 
cells, we need the grand canonical ensemble whose partition function is 



ZG{T,V) = J2^^Z^iT,V) , (2.6) 



e 

N=0 



where /j, is the chemical potential derived from Zjv. The grand canonical partition function 
is related to the equation of state by 

PV 

InZa^ — ^Nil-K) , (2.7) 



where 



1 + 6 



with, as usual. 



W 27iGm^n 



(2.8) 



2K 3T ' ^^(^^^^^ (2-^) 
the ratio of gravitational correlation energy, W ^ to twice the kinetic energy, of peculiar 
velocities in a volume, here taken to be spherical with radius R. The two galaxy correlation 
function is ^2(^)- 

From equations (2.1)-(2.9), the spatial distribution function, which is the sum over all 
energy states in the grand canonical ensemble becomes (Ahmad, Saslaw and Bhat 2002): 
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In the limit e ^ so that C 1 b, f{N, be) becomes the distribution function for 

point masses (Saslaw 2000). 

To derive the gravitational velocity distribution function, f{v), from the spatial dis- 
tribution, f{N,be), we need to modify the grand canonical ensemble by adding a relation 
between N and v. Previously in Saslaw et al (1990), this was postulated by assuming that 
over a given volume any local fliictiiation of kinetic energy (caused by correlations among 
particles) is proportional to its potential energy fluctuation giving 

n( — \ =av\ (2.12) 

\ / Poisson 

The detailed configuration in each fluctuation is represented by the form factor a, indicating 
the local departure from 

1\ ^l£„yrf|U2.18«i (2.13) 



/ Poisson 

for a Poisson distribution (see Saslaw 1987, equation 33.23 with a — 0). Although a varies 
among fluctuations, it clusters strongly around its average value, as n-body simulations show 
(Saslaw et al 1990). 

Instead of postulating (2.12), we can now derive it from the configuration integral in 
(2.1): 



g^(r, V)= f ... f exp[-0(ri, r2, . . . , r^)^- V' 

(2.14) 



= V 



N 
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The local fluctuation of the potential energy around the spatially constant mean field in a 
volume containing N galaxies therefore has an ensemble average value: 

_ 1 dQN _ d{\n Qn) 



Qn ar-i 

= -3(iV-l)6,T ^^-^^^ 
= -{N - l)bem{v^) 

using (2.14) and (2.8). We may also write the average potential energy in this volume as 

W = - E 9^.-911^^/1.) , (2.ia) 
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where 

ji)^ N{N-l) ^^^^ ^ (^)poisson ■ ^^'^^^ 

Here 77 represents the form factor. Prom (2.15)-(2.17), we obtain 

^-^A^V (2-18) 



which is the same as (2.12) averaged over many cells having the same and a with 

26^ 

Gmr] ^ ^ 

In regions of enhanced density, the r^j arc generally smaller than their Poisson values 
for the average density n and consequently r/ > 1 from (2.17). Using natural units G — m — 
R — 1, equation (2.19) shows that for the currently observed clustering with h ^ 0.75, the 
value of a is also of order unity. In these units, averaging (2.18) over the entire system for 
individual galaxies so that (N) — 1 gives the approximate relation with /3 = {v"^) : 



a!v( — ) /3-' . (2.20) 

\ '"ii / Poisson 

Although this average value of a (or equivalently rj) simplifies the description, it gives a 
reasonably accurate result since the velocity distribution involves all the galaxies in the 
system, and is thus an effective average over all the cluster configurations in the entire 
ensemble. 

To transform the spatial distribution f{N) in (2.10) into the velocity distribution func- 
tion, we follow the same procedure as for e = (Saslaw et al 1990) to obtain 

f{v)dv = y^!~y [Q;/3(l - + ahv^r"-^e-^''^^^-''^^+^''^^\dv . (2.21) 
Tyav^ + 1) 

Here P(q;v^ + 1) is the usual gamma function. Pigure 1 shows examples of f{v) for some 
typical parameters suggested by the observations and for a range of e. Por e = 0, we obtain 
the previous results. As e increases, equations (2.3) and (2.8) show that C(e/i?i) decreases 
and that increases monotonically with increasing b and also with increasing ({e/Ri). 
Moreover, b > b^, and the disparity between them will increase as e increases, though for 
e/Ri < 0.5, the value of b — b^ < 0.1. As e increases. Pig. 1 shows that the peak of f{v) shifts 
towards higher velocities, but there are fewer galaxies in both the low and high velocity tails, 
relative to the point mass case e = 0. This effect only becomes substantial for e > 0.5, when 
the haloes of individual galaxies overlap with the haloes of their nearest neighbours. 
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The effects of halo overlap are to smooth out the gravitational potential wells so that 
close encounters produce fewer high velocity galaxies. The reason that there are also fewer 
low velocity galaxies as e increases is that the forces are more uniformly distributed through- 
out the system and fewer galaxies are so gravitationally isolated that their peculiar velocities 
simply decay adiabatically as the universe expands. 

These effects are also illustrated nicely by the radial velocity distribution function, 
which is what we usually observe. The total velocity v is related to its radial and tangential 
components by 

= vl + vl ^vl = vl + v\ . (2.22) 

Integrating (2.21) over the tangential velocities (Inagaki et al 1992) now gives the radial 
velocity distribution function for the softened potential: 

j{vr) = o?i5{\ - 6,)e-"^(^-''^) 

^ r v± W - be) + aKiv'r + ^1)]"^"-+"^^-^ ) , (2-23) 

Jo ^MT^ r[a{v^^ + vl) + 1]] ^ 

where again /3 = {v^ + v\). 

Fig. 2 shows the observed histogram of radial peculiar velocities for a sample of 825 Sb- 
Sd galaxies with t^radiai ^ 5000 km s~^ from Figure 2a in Raychaudhury and Saslaw (1996), 
as an illustration. The solid line is the best fit to this histogram for e = 0; it has a = 40.9, 
f3 = 1.54 ^ 3(f,^), and b = 0.91. The broken lines show the effect of increasing e/Ri while 
holding the other parameters constant. 

Evidently only relatively large haloes, e/Ri > 0.3, affect the velocity distribution (or 
the counts in cells) significantly. Since these large haloes soften and reduce the gravitational 
forces, they lead to reduced clustering. This implies fewer galaxies in the very high (radial) 
velocity tail of the distribution and a broadening around the most probable velocity (at zero 
for Vr) whose peak is shifted to somewhat higher velocities. This leaves fewer underdense 
regions with low velocities. If there are fewer rich clusters and the number of galaxies is at 
least approximately conserved (i.e. mergers do not dominate, of if they do , they mostly 
occur near the centers of mass of the merging galaxies and these centers have approximately 
the same distribution as the galaxies themselves) then fewer voids and underdense regions 
are formed to create the clusters. 

As an illustration, if we take Ri to be an average distance of ~ 1 Mpc between typical 
massive galaxies, this would suggest that haloes would need to be about 300 kpc in radius 
for a significant effect. These numbers are made uncertain by the criterion adopted for 
massive galaxies and the lack of precision in determining the "effective radius" of a halo 
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whose density falls off as r~^, but they are indicative. If e/Ri ^ 0.5, the haloes are nearly 
touching, and if e/i?i > 1, several galaxies would share a common halo. The latter case, 
although it is favoured by some cold dark matter scenarios, would not appear to agree with 
the presently observed velocity distribution. Instead of the orbits being mainly determined 
by the (softened) gravitational forces of individual galaxies, such models with e/Ri > 1 
determine the orbits mainly by the inhomogeneous distribution of the matter in the common 
halo. For a communal halo which was not highly concentrated, this would lead to a radial 
velocity distribution which is relatively broader and less peaked compared to the case with 
individual haloes and small e/Ri. 

In Fig. 2, by setting e/Ri = 0, we determine the best fit values of b, a and p. Including 
effects of dark matter haloes, observations of the spatial distribution function, equation 
(2.10), or the variance of its counts in cells, gives be directly. Observations of ^2('^), the 
peculiar velocity dispersion and an estimated average galaxy mass m (including haloes) would 
provide b from equation (2.9). Then equations (2.17) and (2.20) give a and all the quantities 
in the theory can, in principle, be found consistently from observations. Comparing be and 
b yields e/Ri. 

Figure 2 in Ahmad, Saslaw and Bhat (2002) shows that for e/Ri < 0.5, the values of b 

and be differ by less than 0.07. Although this difference is small, our Figure 2 shows it can 
significantly affect the peak of the radial velocity distribution. In practice there is not yet 
enough information to determine e/Ri accurately, and it can only be fitted within ranges 
for prescribed confidence limits. Current data, consistent with e/Ri < 0.5, suggest that 
most individual galaxies are surrounded by their own haloes rather than embedded in large 
massive common haloes. 



3. PROBABILITIES FOR BOUND AND VIRIALIZED CLUSTERS 

3.1. The Density of Energy States 

For estimates of the mass-luminosity ratio or the amount of dark matter in clusters of 
galaxies, it is generally assumed that the clusters are gravitationally bound and have a high 
probability of being virialized (Saslaw 1987, 2000). Our statistical mechanical analysis now 
makes it possible to calculate these probabilities directly for the cosmological many-body 
problem, and to determine how they are affected by the value of b and the average peculiar 
velocity dispersion (or temperature), as well as by the average halo size e/Ri. 

The general probability density, in a grand canonical ensemble, for a region with N 
galaxies to be in the differential range of energy states, dE, follows from the first equahty of 
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equation (2.10) without the summation over energies as 

P{E, N, To, V, i,)dE = g{E)^-^^;-y-^dE (3.1) 

where g{E) is the density of energy states and Zq is given by equation (2.7). The continuous 
density of states arises from transforming the sum over individual energy microstates Ei 
into an integral over a finite range of energies. Here To is the average temperature for the 
ensemble. 

To determine g{E) for cells containing N galaxies, we first recall, for example from 
Greiner et al (1994), that a subset of the grand canonical ensemble in which each cell contains 
exactly N galaxies (particles) forms a canonical ensemble. It is in contact with a temperature 
reservoir so that energy, but not particles, can move among cells. A particular cell of volume 
V , with a particular configuration of N galaxies has an energy E. Averaging these energies 
over all cells gives the average energy U which appears in the equation of state for the system 
(Saslaw and Hamilton 1984; Ahmad, Saslaw and Bhat 2002) 

U = ^-NT{1 - 2b,) , (3.2) 

since if all cells have a given value of A^, the grand canonical equation of state also applies to 
this canonical ensemble. Boltzmann's fundamental postulate of statistical mechanics relates 
g{E) to the entropy in the microcanonical ensemble by 

n{E)^ f g{E)dE ^ e^^""'^'^^ (3.3) 
Jae 

where Q{E) is the number of energy states in a small range AE <^ E, so Q(E) = g{E)AE. 

It is well known that for large values of A^, the energy E in the microcanonical ensemble 
is very nearly equal to the average energy U in the canonical ensemble. Under this condition, 
we can examine canonical ensembles with different average energies U to determine P{E) 
from equation (3.1). Therefore we next determine the detailed condition for S{E,V,N) 
S{U, V, N) to hold. 

The energy dependence in equation (3.1), which involves both the density of states and 
the Boltzmann factor, is the free energy —F — T S{E) — E. The probability of any energy 
will therefore have an extremum at 

rm ^i. (3.4) 



dEJ 



E=U 



This extremum will be a maximum if 



d''S \ _ fd_l\ _ _J_ /dT\ _ 1 
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where Cy is the specific heat of the entire ensemble. This ensemble will contain a range 
of cluster masses and energies for any particular value of b^. As equation (3.7) indicates, 
Cy will remain positive for less than 0.86. If Cy becomes negative, P{E) will have no 
maximum. This state represents a system dominated by negative energy clusters which are 
however virialized so that they do not become unstable on short timescales and the system 
remains in quasi-equilibrium. When there is a maximally probable state for energy E, we 
expand the free energy around it: 



S{E) - § 



S{U) - ^ 
S{U) - \ 



V 



(3.6) 



Hence the departures of S{E) and E from S{U) and U are closely related to the variance 
of the energy fluctuations around the mean. Using earlier results for Cy and {(AE)"^) in 
the grand canonical ensemble (Saslaw and Hamilton 1984; Ahmad, Saslaw and Bhat 2002), 
which give an upper limit to the fluctuation term in equation (3.6), we flnd 

1 iry _ 1 (5 - 20^. + 34&g - 166^) _ 



2NT^Cv^ ' 4(l-6,)2 (1 + 46,-662) 



This ratio has two poles at 6, = (2 + vl0)/6 = 0.86 when Cy passes through zero and at 
b^ — \ when all the galaxies collect into one nearly virialized cluster (Baumann, Leong and 
Saslaw 2003). Away from these poles, the fluctuation term decreases rapidly (e.g. it is 44.5 
at 6e = 0.8, 15.2 at b^ = 0.75 and 3.78 at bf = 0.5) compared to S{U) — U/T which is of 
order A^[l + 36 + ln(l — b)+ non-gravitational terms]. Since groups must contain N > 2 
galaxies, the approximation S{E) ^ S{U) is usually reasonable. This is closely related to 
the neglect of terms of order In N in the thermodynamic limit of large N since equation (3.3) 
then gives S — In g{E) + In A£' In g{E). Taking AE to be a unit dimensionless variation 
of the energy shell, AE^, from equation (3.10) below gives S — In g{E^). 

The grand canonical entropy, which for a sub-ensemble of fixed is the same as the 
canonical entropy, is a function of temperature (Saslaw and Hamilton 1984; Ahmad, Saslaw 
and Bhat 2002): 

SiT,V..,.-m.{JL^),mnl,^ aT-, _ , 5, , 2^ ,„ (^) (3.8) 

where from equations (2.2) and (2.9) 

a^liGmrncd-) (3-9) 



2' ' V^i 



- 12 - 



has dimensions of (energy)^. In fact a}/^ is essentially the absolute value of the potential 
energy of two galaxies with average seperation. The relation between T and E for U ^ E 
in volumes with energy E and given is given by the energy equation of state (3.2). 
Normalizing to dimensionless variables, 

277 

E, = r , (3.10) 

3A^a3 

and 

= — . (3.11) 

as 

Here is defined as the local fluctuating temperature within the ensemble, which differs 
from To* = 7o/a^/^ the average dimensionless temperature of the whole ensemble. Using 
equation (2.8), equation (3.2) becomes 

E* - ' (3-12) 

or 

- E,T^ -T,-E, = 0. (3.13) 

The solutions of both the linear equation for E^{T^) and the quartic for T^^^E^), although 
equivalent, provide different insights into the relation between T and E. 

Since equation (3.12) is simpler, we consider it first. Figure 3 shows E^[T^). It is double 
valued between the zeroes at = and = 1 and has a minimum at 

T, = - 1^' = 0.54, E^ = -0.390 . (3.14) 

This value of will be found to divide the two real branches of the solutions to the quartic 
equation (3.13). For large values of T*, produced either by a high temperature or a weak 
interaction a^/^ 0, E^ = in equation (3.13), which becomes the perfect gas equation of 
state. For small T^, the limit of equation (3.13) is E^ = — which is the equation of state 
for a completely virialized system having specific heat Cy — —3/2. 

The real solutions of the quartic equation (3.13) are 



E 1 / 1 
n^{E^) ^^ + -^^ + Z,±-^ (3.15) 



where 




, ^ /I 223£;3 Ef\ - 1 El I 223EI E^ . 
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and 

The number of energy states from equations (3.3) and (3.8) is 

) ' (^) + ^^''^"^ e-^^/(^+^- ) . (3.18) 

We then find T^,{E^) from either equations (3.12) or (3.15)-(3.17) to calculate g{E^) = 
dQ{E^)/dE^ ^ Q{E^:)/AE^:. In fact, we use the results from both these solutions to double 
check the numerical computations; they agree with each other. Already from equation (3.18), 
we see for > 1 and thus ^ E^ that g(E.,) oc ^^^^ ^ which agrees with the usual 
energy dependence of g{E) in an ideal gas (Greiner et al 1994). For the dual limit as T,, — » 
in equation (3.13) and — > —E^, recalling that the potential energy dominates so that 
< 0, we obtain the same energy dependence of g{E^). This describes a dense or massive 
system in which a^^^ oo, even though T may remain large and virialized, and the system 
behaves approximately as an idealized gas gravitationally bound in a cluster. A system of 
approximately isothermal spheres each with the same nearly Maxwell-Boltzmann velocity 
distribution would be an example. In equation (3.18), we may normalize the volume V and 
the mass m to unity by selecting the phase space cell normalization A^, as often done for 
the classical gas (Sommerfeld 1955). 

Generally we solve for g{E^) numerically from (3.18) and either (3.12) or (3.13) and then 
integrate equation (3.1) over a range of energies (or equivalent temperatures) of interest. The 
normalization is such that for fixed N , Tq*, V and fi 



f 



-^*max 

P(E„ iV, To,, V, ij)dE, = 1 . (3.19) 



It is important to keep in mind, because T^{E^) has a double valued regime as in Figure 
3, that the subscripts "max" and "min" refer to the maximum and minimum values of 
over which the integral is evaluated, and not to the maximum and minimum values of the 
energy. Thus the integral is evaluated from a value of E^. min = E^{T^ min) corresponding to 
the minimum value on the branch, through E^ viriai = —0.33 on the branch, and 
then to £'*max = E^{T^ max) after switching to the branch at = 0.54 and E^ — —0.390. 

Since the system is in quasi-equilibrium in the expanding universe, its range of energies, 
-E'*min <-£'*< ^'^max, is morc restricted than the usual infinite range. States with positive 
energies represent unbound groups or clusters, and if their energy is too great they will 
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fly apart on time-scales much shorter than the Hubble time and be inaccessible to the 
thermodynamics of the quasi-equilibrium system. Negative energy states which are not 
near virial equilibrium will collapse on time-scales much shorter than the Hubble time, 
and they will not be accessible either. We can estimate -E^viriai directly from £'*viriai ~ 
-E'viriai/(l-5A^a^/^) ~ i?viriai / 1 | ~ —0.33 ou the branch of the two real solutions of the 
quartic equation. One can also obtain a less direct but similar estimate of £^*viriai from the 
temperature r^,vii.iai of a cluster modeled as a polytrope. To estimate the maximum kinetic 
energy, we note that if the kinetic energy is twice the potential energy, the velocities in the 
unbound region will be about twice the virial velocities that the region would have had if it 
were viriahzed. These would be about the escape velocity, giving E'^max ~ 1- More precise 
values of these limits will follow from the solutions for (?(£'*) and the binding probabilities. 
In Figure 3, we show the regions of these different solutions. Two of the four solutions of the 
quartic equation become complex for < —0.39 which is where the switch from the to 
the T*_ real branches occurs. For > 1, the total energy is positive and tends towards the 
perfect gas limit. With these solutions for T^{E^), we obtain the density of states 



dE^ dT^ dE^ 

^^-N\nN+5N/2 (^l J^-^^N ^-3N/{1+TS) j^{3N/2)-1 



(3.20) 



represented here in terms of for simplicity, and setting 

= 1 . (3.21) 



Prom (3.20), we immediately see that there is a smooth transition around — where the 
total energy becomes negative. As — > 0, the density of states becomes proportional to 
T* ^-'^Z^. This infinite number of negative energy states describes a singularity. To reach this 
singularity, the motions must be dominated by radial infall with negligible random velocities 
and therefore negligible temperature. In this regime, the quasi-equilibrium nature of the 
theory breaks down. 

Figures 4 and 5 illustrate g{T^,) as a function of T*. which is single valued unlike E^, and 
also for a range of N. It is clear from the linear plots in Figure 4, and from the logarithmic 
plots in Figure 5 for a much greater range of N, that the (?[T* (£"*)] have infinite maxima 
both for the perfect gas and for gravitating systems of galaxies as T^, — > 0. In the latter limit, 
the random velocity dispersion {v"^) 0, or the gravitational energy ~ a^/^ — > oo, and the 
galaxies in a cell collapse into a singularity having infinite entropy. The dynamical timescale 
for such a collapse is so short that these systems cannot be in quasi-equilibrium, and such 
states are inaccessible to the statistical thermodynamic description. Moreover, they are also 
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inaccessible because they are not physically realistic. In any collapse toward a singularity, it 
is very unlikely that all the radially infalling orbits will remain in phase and reach the centre 
simultaneously. Instead, they perturb each other, introducing more random velocities which 
build up until the contracted system reaches virial quasi-equilibrium which then evolves on 
a much longer timescale. As shown in §3.2, this occurs around T* ^ 0.1. Between this value 
and T* 1.3 the density of energy states has a broad minimum corresponding to states 
which have neither a high gravitational entropy nor a high perfect gas entropy. The total 
number of configurations of these states is constrained by their relatively narrow range of 
energy A£'* for a given range of temperature AT* around the minimum of the E^{T^) relation 
in Figure 3. Therefore these are much less probable low entropy states. 

Figures 4 and 5 also show the somewhat anomalous behaviour of the case N — 2, 
an effect of the approximation 3N/2 {3N/2) — 1. However, for > 4, the curves 
become more regular and this approximation becomes more reasonable; even for = 2, it 
is qualitatively accurate. For larger N. the transition to large g becomes sharper and occurs 
at both lower and higher values of T* than for small A^. With g{E^), we can now determine 
P{E^) for a canonical ensemble with each cell containing A^ galaxies using equation (3.1). 
The normalization (3.19) cancels the contributions of A^ in the chemical potential and Zq 
although N still affects the form of g[E^,{T^)\. 



3.2. Binding and Virialization 

The probability that A^ galaxies in a cell are gravitationally bound together is 

Phound{N,To^) = ^r.'""" ^ . (3.22) 

J^'-''P{E,,N,To.)dE, 

The number of bound states determines the number of unbound states through the normal- 
ization Pbound + -^unbound — 1- We havc also computcd i-*unbound Separately as a check that 
the integrations satisfy this normalization. 

As figures 4 and 5 show, the integrals depend somewhat sensitively on the limits £^*min 
and £'*max through their dependence on Tlmiin and Tl^max ■ We will explore this sensitivity 
around average values of these limits. To estimate an average value for E^^max, we note that 
for the i-ih cell containing A^ galaxies 







Wi\ 




\W 





E,i = ' - , " (3.23) 



where = a^/^ is the average gravitational correlation energy. If |Wj| = in a typical 
state, then if the cell were bound in virial equilibrium it would have Ki — \W\/2. If it were 
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just bound, Ki — \W\, and if it were sufficiently unbound that the root mean square velocity 
were twice the virial root mean square velocity, so that Ki = 4ii'viriai = 2|Ty|, then the 
configuration would expand to about twice its size on a dynamical time scale. Such states 
are not in quasi-equilibrium and would not be accessible to the statistical mechanics or 
thermodynamics of the ensemble. Therefore wc exclude them. (In the statistical mechanics 
of perfect gases, this exclusion is not necessary since particles of all energies are assumed to 
be in equilibrium.) This suggests -E^max ~ 1 on the T*_|_ branch. 

To estimate i?*mirn it is more straightforward to work with T*min, which is a positive 
single-valued quantity, and then find i?*min through equation (3.12). Using equations (3.9) 
and (3.11) and the relation between temperature and kinetic energy, gives 

^_ n'^Hv') _ 0.0881 
* 3(3/2)1/3(1/3^/5 (1/3 (ri)2 ■ ' 

Here we have used the relation (cf equation 33.22 in Saslaw 1987) connecting the average 
uniform number density n to the average seperation (ri) between galaxies: n^/^ = 0.55/ (ri). 
The dynamical timescale is r = (Gp)~^/^, and (^/^ is usually unity or slightly less from 
equation (2.3). Although the minimum value of the velocity dispersion for quasi-equilibrium 
to apply is not precisely defined, it is necessary that (t'^)r^/ (ri)^ ^ 1 so that the configuration 
does not collapse on a dynamical time scale. This suggests T^min ~ 0.1 corresponding to 
-E*min ~ —0.1 on the T*_ branch. 

To check these estimates, we can compute Phound{N, E^j^i^) from equation (3.22) for a 
range of -E*min- The binding probability will also depend on through its relation to Tq^, in 
equation (2.8). Figure 6 shows the results for —0.12 < E^^min < —0.04 and 6^ = 0.75. There 
is a bifurcation at £'*min = —0.10 between bound probabilities which increase with N and 
those which decrease. We would expect that the more galaxies there are in a cell of given 
volume, the more likely they are to form a gravitationally bound group. Therefore functions 
-Pbound(A^) which dccrcasc with increasing are unphysical and, in Figure 6, only those cases 
with increasing Pbound(A^) are valid. These unphysical cases have -E*min < —0.10. However, 
E*min ~ —0.1 is also the smallest value from equation (3.24) which is consistent with quasi- 
equilibrium. Both arguments therefore suggest that -E*min ~ —0.10 on the T*_ branch of 
the E^{T^) relation in Figure 3, is the physical lower limit of the probability integrals for 

-fbound(-^)- 

With E^-cain = —0.1 on the T^,_ branch. Figure 7 shows the effect of on the binding 
probabilities for different A^. Note that N enters through a^/^ in equations (3.9)-(3.11) and 
the normahzation in equation (3.21). As increases from 0.1 for a nearly perfect gas to 0.9 for 
a highly clustered gravitational system, the binding probabilities increase significantly. These 
probabilities are high even for relatively small values of 6^ which illustrates the efficiency of 
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gravitational clustering. Dynamically, positive density perturbations expand more slowly 
than the surrounding universe even if they are linear, and eventually they become bound. 
In a statistical mechanical description, these bound states have many more possible velocity 
and density configurations than the corresponding unbound states. 

One range of bound energy states is of special interest. These are the approximately 
virialized states around ^ —0.33 on the T^,+ branch, discussed after equation (3.19), 
which is 2f/^,/3 in equation (3.10) with U^, = —1/2. Although the definition of the "virial 
range" is somewhat arbitrary, here we take it to be —0.39 < < —0.32 on the T^j^ branch. 
Recall that for E^ < —0.39, the T^{E^) relation becomes imaginary. Figure 8 shows bound 
and virialized probabilities as well as the conditional probability that a bound system of 
particles is in this virial range. Most bound states become virialized and then evolve slowly 
so we should expect that many, including our local group have been observed. 

3.3. Compcirison with Observed Groups of Galaxies 

Identifying physically bound groups of observed galaxies is often uncertain based on 
limited velocity and positional information. Nevertheless, reasonable estimates of the num- 
bers of such groups have been made for the local universe where problems of incomplete 
sampling have been minimized. These estimates also depend on the criteria used to define 
a group. Garcia (1993) has explored these issues in some detail and developed a catalog of 
485 groups with > 3 members each in a selected sample of 6392 nearby galaxies with 
apparent magnitudes Bq < 14.0 and radial velocities < 5500 km s~^. Two procedures, one 
based on hierarchial structures, the other based on percolation, are used to establish group 
membership. From these, Garcia (1993) derived the number of groups having A^ > 3 mem- 
bers shown in Figure 9 as the solid (hierarchial) and dashed (percolation) histograms taken 
from Garcia's Figure 2. The moderate difference between these histograms is a measure of 
the uncertainty of group membership. We use these results for comparison with the expected 
probabilities of bound groups in cosmological gravitational many-body clustering. 

To calculate the theoretical numbers of groups at the present time (using the ob- 
served 6e = 0.75 so that To* = 3~^^^ — 0.70), we multiply f{N) from equation (2.10) by 
Pbound(-^) 2o*) from equation (3.22) and normalize the total number to the 485 groups in 
Garcia's catalog. Figure 9 shows the results with -E^min = —0.10 for the solid line. We have 
also examined the results of varying -E'*min between -0.04 and -0.1. This range of £'*min does 
not affect the probabilities significantly. 

Overall there is reasonable agreement between the theoretical estimate of the number 
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of bound groups and the observed estimate for 3 < A?" < 10, which includes most of the 
sample. For 10 < iV < 17 the predicted number is higher, while for = 18, 19, it agrees 
again with the observations. The causes of the discrepancy are difficult to determine, but 
may be related to the sensitive dependence of the numbers in more populated groups on the 
magnitude cutoff. The theory assumes that all galaxies have the same mass and luminosity 
for simplicity and this could be improved. Larger catalogues, complete to fainter limiting 
apparent magnitudes will soon become available and it will be interesting to see if they 
improve the agreement. 



4. THE COOL LOCAL HUBBLE FLOW 

We can apply these results to the long-standing problem of the small peculiar velocity 
dispersion, < 60 km s~^, of galaxies in our local group (van den Bergh 1999, 2000; Ekholm et 
al 2001). Consequently, the relatively small mass of the local group, (~ 2±0.5 x IO^^Mq) does 
not significantly perturb the linear Hubble flow beyond about 1.5 Mpc (Sandage, Tammann 
and Hardy 1972; Sandage 1986; Ekholm et al 2001). As Sandage (1999) has emphasized 
"... the explanation of why the local expansion fleld is so noiseless remains a mystery." 
The dynamics of the local flow and more general catalogs of peculiar velocities have been 
explored by many authors (Groth, Juszkiewicz and Ostriker 1989; Peebles 1990; Burstein 
1991; Governato et al 1997; Willick et al 1997; Wilhck and Strauss 1998; van der Weygaert 
and Hoffman 2000; Nagamine et al 2001). 

Previously, there have been three main types of explanation for the cool local Hubble 
flow. One, pioneered by Kahn and Woltjer (1959) and subsequently often discussed, involves 
specific detailed models of the dynamics of the formation of the local group. On the other 
hand, there may be more generic reasons for low peculiar velocities. Sandage, Tammann 
and Hardy (1972) suggested two possibihties: a very low density universe, or a high density 
universe dominated by a uniform component of matter. In both these cases, gravitational 
galaxy clustering would usually be relatively weak if it started from weakly clustered initial 
conditions. The low density universe is not generally considered likely at present, but the 
uniform high density possibility has recently been revived by Barychev et al (2001) in terms 
of a dark matter component suggested by a large positive value of the cosmological constant. 

In explorations of these or other explanations, the fundamental question is whether our 
weakly clustered local group is unusual enough to demand a specific dynamical history, or 
whether it is quite commonplace. Therefore we need to calculate the probability that a small 
weakly bound group such as our local group can form. Although our local group contains 
about three dozen galaxies in a volume of about one megaparsec radius, most of its mass 
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and energy are dominated by Andromeda and the Milky way (including their dark matter 
haloes). Most of the other galaxies are dwarfs or satellites. This makes us a group with 
effectively iV = 2, or perhaps N — 3 major galaxies if all the smaller galaxies are equivalent 
to one large galaxy. 

Prom Figure 6, at least 86% of all such groups are bound. Thus the probabihty of a 

bound group such as our local cluster is approximately > 0.86 f{N,b^). Using equation 
(2.10) for f{N, h^) with = 1 and = 0.75 gives lower limits to the probability that there 
are bound groups with N effective massive galaxies within a radius of one megaparsec from 
an arbitrary point in space. These limits are 0.033 and 0.019 using effective values of = 2 
or A^ = 3 respectively for our local cluster. Thus within a radius of 10 Mpc, we would 
expect more than about (47r/3) (1000) (0.033) = 135 for A^ = 2 or 80 for A^ = 3 groups such 
as our own whose masses are too small to perturb their local Hubble flow strongly. Prom 
equation (2.10), there would also be bound groups of greater mass which would perturb the 
Hubble flow more strongly, but such massive groups would be less numerous. Low mass 
bound groups like ours should be separated on average by ~ 10(80)-^/=^ ^ 2.3 Mpc. So if 
they have a radius of ~ 1 Mpc and are randomly placed there is about a 10% probability 
that a galaxy at a random position of space will belong to such a group. Many galaxies will 
belong to even lower mass groups. This is also a lower limit since most groups and galaxies 
are clustered rather than located randomly. Therefore there is at least a 10% chance that an 
astronomer in an arbitrary galaxy would belong to a group which does not seriously disturb 
the surrounding cosmic expansion. This seems to us a high enough probability to remove 
most of the mystery of the cool local Hubble flow. 



5. SUMMARY AND DISCUSSION 

Distribution functions for peculiar velocities and gravitational binding probabilities 
among galaxies provide new insights into their clustering. These complement the more 
often analyzed distributions and correlations of galaxy positions. Galaxy haloes, for exam- 
ple, generally affect velocities more strongly than they affect positions. Haloes soften the 
gravitational forces and inhibit the high velocities associated with strong clustering. Strong 
clustering can still develop at these lower velocities, it just takes longer. 

Gravitational statistical mechanics is proving to be a powerful method for calculating 
these distributions. It applies not only to the cosmological many-body problem, but also 
to finite systems (de Vega and Sanchez 2003), although since finite systems have different 
symmetries they have a different thermodynamic limit than our infinite system. In our 
cosmological case, the statistical mechanical derivation of the peculiar velocity distribution 
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function does not require assuming that average local kinetic energy fluctuations are pro- 
portional to average local potential energy fluctuations. We show that this assumption is a 
direct consequence of the partition function. Moreover our new derivation generalizes pre- 
vious results by including possible haloes around galaxies. The currently observed peculiar 
velocity distribution function appears most consistent with each massive galaxy usually be- 
ing surrounded by a single individual halo rather than many galaxies sharing a large common 
halo. Future extensive and accurate radial peculiar velocity observations can quantify this 
result more precisely. 

The density of energy states for a system of gravitationally interacting galaxies in the 
expanding universe is a fundamental quantity. Usually the classical density of energy states 
can be calculated in detail only for non-interacting systens, but our cosmological gravitating 
system turns out to be a rare exception. Our calculation of this density of states makes it 
possible to flnd the probability that a group oi N >2 galaxies is gravitationally bound and 
virialized. This calculated probability agrees well with the observed probability distribution 
for flnding groups of N galaxies within about 100 Mpc. The preponderance of relatively low 
mass groups also provides a simple explanation for the observation that most groups, such 
as our own, do not disturb the Hubble flow beyond the group appreciably. 

These results also show that the efficiency of gravitational galaxy clustering is very high. 
More than 85% of double systems and more than 95% of groups with > 10 massive galaxies 
arc bound. Of these at least 95% arc virialized. This is consistent with earlier studies (Saslaw 
1979, 1987) which found that the observed form of the galaxy two-point correlation function 
leads to highly efficient clustering. 

The binding probability increases rapidly as N increases. This is more dramatic for 
small values of 6^ than for large values. For large b^, nearly all groups are bound for N >2. 
Since increases as the universe expands (a manifestation of entropy increase), this describes 
the evolving formation of bound groups and clusters. At the present time, nearly all clusters 
with > 60 massive galaxies are expected to be bound and virialized. Even for moderate 
values of 0.25, which occur at redshifts ~ 2.5 — 5 depending on the details of the 

expansion rate (Saslaw and Edgar 2001), we would expect most groups to be bound and 
virialized. 

The high efficiency of binding and virialization on relatively small scales provides further 
insight into the quasi-equilibrium nature of galaxy clustering. On small scales, virialized 
clusters remain unstable only over dynamical relaxation timescales ~ N/^/Gp, which are 
generally longer than a Hubble time. On large scales, the formation of new clusters from 
linear perturbations takes at least a Hubble time. Therefore changes in b and the global 
equations of state take at least a Hubble time for their instability. This exceeds the local 
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dynamical crossing timescale, l/\/Gp, and is thus consistent with the approximation for 
quasi-equihbirum (Saslaw 2000). 

Since the probabihty that a group forms depends on its total energy, and its energy 
depends on its shape as well as on its size and number of massive galaxies, our results make 
it possible to calculate the probability that stable or unstable structures such as filaments 
or "great walls" can occur. We will analyze this elsewhere. 
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the cool local Hubble flow and groups of nearby galaxies. We also thank Phil Chan, Alvaro 
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Fig. 1. — Peculiar velocity distribution functions using observed values of a, /? and b (Ray- 
cliaudliury and Saslaw 1996) from the sample of 825 Sb-Sd galaxies with distance D < 5000 
km s~^ from the Mathewson et al catalogue with peculiar velocities corrected for the Dressier 
et al. (1986) bulk motion of 599 km s~^. The velocity distributions of equation (2.21) are 
compared for a range of values of the halo size, e/ Ri. 
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Fig. 2. — The observed distribution (histogram, Raychaudhury and Saslaw 1996) of radial 
pecuhar velocities for the sample of 825 Sb-Sd galaxies with distance D < 5000 km s^^ 
from the Mathewson et al catalogue, with peculiar velocities corrected for the Dressier et al. 
(1986) bulk motion of 599 km s~^. The radial velocity distributions of equation (2.23) are 
compared for a range of values of the halo size, e/Ri. 
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Fig. 3. — The relation between and The curve corresponds to the two real branches of 
the four solutions in the quartic equation (3.13). The curve reaches a minimum at T^, = 0.54 
and = —0.39. The region < 1.0 represents bound groups. The solid line represents 
the T*_|_ branch and the dashed line represents the T*_ branch having real valued energy. 
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Fig. 5. — The logarithm of the density of energy states for various values of A^. The curves 
all reach a minimum point at T*=0.54. 
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Fig. 6. — The bound probabilities for —0.12 < E^rain < —0.04. 
equilibrium limit. 



—0.1 is the quasi- 
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Fig. 7. — The effect of 6^ on tlie bound probability. As increases, tlie bound probability 
of clusters with N >2 increases, indicating greater efficiency in gravitational binding. 
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Fig. 8. — The probabilities of bound and virialized clusters with N > 2 galaxies. The 
conditional probability that a bound cluster is virialized is also plotted. Values of — 0.75 
and E^rnin = — 0.1 are used. 
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Fig. 9. — Histograms are from Figure 2 in Garcia (1993) showing the observed number 
distribution of nearby groups using the hierarchial method (sohd hne) and the percolation 
method (dashed hne). The solid curve with ftg = 0.75, A'^ = 1.0 and -E*min = —0.1 is our 
theoretical result. 



